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f"**. ■ Abstract 
(N " 

Starting from the spectral analysis of g-circulant matrices, we consider a new multigrid 
method for circulant and Toeplitz matrices with given generating function. We assume that the 
^ size n of the coefficient matrix is divisible by g > 2 such that at the lower level the system is 

, reduced to one of size n/g by employing g-circulant based projectors. We perform a rigorous 

two-grid convergence analysis in the circulant case and we extend experimentally the results 
to the Toeplitz setting, by employing structure preserving projectors. The optimality of the 
proposed two-grid method and of the multigrid method is proved, when the number 6 £ N of 
recursive calls is such that I < 6 < g. The previous analysis is used also to overcome some 
pathological cases, in which the generating function has zeros located at "mirror points" and 
the standard two-grid method with g = 2 is not optimal. The numerical experiments show the 
' correctness and applicability of the proposed ideas both for circulant and Toeplitz matrices. 

o 
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1 Introduction 

In the last 20 years multigrid methods have gained a remarkable reputation as fast solvers for struc- 
tured matrices associated to shift invariant operators where the size n is large and the system shows 
a conditioning growing polynomially with n (see [TDl EH EH El El E51 El [H] and the references 
therein). Under suitable mild assumptions, the considered techniques are optimal showing linear 
or almost linear (0(n log n) arithmetic operations as the celebrated fast Fourier transform (FFT)) 
complexity for reaching the solution within a preassigned accuracy and a convergence rate indepen- 
dent of the size n of the involved system. These excellent features carry over the multilevel setting 
and mimic very well those already known in the context of elliptic ordinary and partial differential 
equations (see [HI [THl Ell [22] and the references therein). In particular, if the underlying structures 
are also sparse as in the multilevel banded case, then the cost of solving the involved linear sys- 
tem is proportional to the order of the coefficient matrix with constant depending linearly on the 
bandwidths at each level. We mention that the cost of direct methods is 0(n log n) operations in 
the case of trigonometric matrix algebras (circulant, r, . . . ) and it is 0(n <* ) for d-level Toeplitz 
matrices (see |16j). Concerning multilevel Toeplitz structures, superfast methods represents a good 
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alternative, even if the algorithmic error has to be controlled, with a cost of 0(n~3~ log 2 (n)): the 
cost is really competitive for d = 1, while the deterioration is evident for d > 1 since it is nontrivial 
to exploit the structure at the inner levels (see [7J and references therein). Moreover, in the last 
case the most popular preconditioning strategies by matrix algebra can be far from being optimal 
in the multidimensional case (see [21] )■ On the other hand, multigrid method are optimal also for 
polynomially ill-conditioned multidimensional problems and they can be extended to the case of 
low rank corrections of the consider structured matrices, allowing to deal also with the modified 
Strang preconditioner widely used in the literature (see [6] and the references therein). 

The main novelty contained in the works from the structured matrices literature is the use 
of the symbol. Indeed, starting from the initial proposal in |10| . we know that the convergence 
analysis of the two-grid and V-cycle can be handled in a compact and elegant manner by studying 
few analytical properties of the symbol (so the study does not involve the entries of the matrix 
and, more importantly, the size n of the system). Already in the two-grid convergence analysis, it 
is evident that the optimality can be reached only if the symbol / has a finite number of zeros of 
finite order and not located at mirror points: more explicitly, if xq is a zero of / then /(xq + ir) 
must be greater than zero. Here we show that the second requirement is not essential since it 
depends on the choice of projecting the original matrix of size re into a new one of size re/2. The 
latter is not compulsory so that, by choosing a different size reduction from n to n/g and g > 2, 
we can overcome the pathology induced by the mirror points. Other approaches for dealing with 
such pathologies were proposed in [5j ITS]. 

In this paper we propose a new multigrid method where the fine problem of size n is projected 
to a coarser problem of size n/g, g = 2,3, ... . We perform a two-grid analysis using the ideas in 
|23j for circulant structures and by exploiting the spectral analysis of g-circulant matrices already 
performed in [18] . As shown in [9], such two-grid analysis is an algebraic generalization of the 
classical local Fourier analysis and it allows one to apply the results obtained for circulant matrices 
also to Toeplitz matrices. A feature of our multigrid is that the coarse problem of size n/g with 
g > 2 allows to obtain optimal multigrid methods with g — 1 recursive calls: it is enough to perform 
the analysis of the arithmetic computational cost related to the size reduction n/g between two 
consecutive levels and to invoke the results in [27]. A further property of the proposed multigrid 
is that the pathologies induced by the mirror points are bypassed as previously described. Our 
proposal could be extended to the multilevel case by tensor product arguments considering the 
increasing of the number of "mirror points". Moreover a V-cycle convergence analysis could be 
performed by following the steps in [21 [1] as a model. A rigorous study in this directions will be 
the subject of future research. 

The paper is organized as follows. In Section[2]we set the problem by recalling the main features 
of the circulant, Toeplitz and g-circulant matrices. In Section [3] we report definitions and classical 
convergence results concerning two-grid and multigrid iterations from |19|. [27] . In Section |4] we 
define our grid-transfer operators and we study the properties of the coarse matrix obtained by 
the Galerkin approach. Section [5] is devoted to the proof of convergence of our multigrid when 
applied to circulant matrices and to briefly discuss the pathologies, that are eliminated by our 
algorithmic proposal. Section [6] is concerned with numerical experiments regarding circulant and 
Toeplitz matrices (e.g. ill-conditioned linear systems coming from approximated differential and 
integral problems). Section [7J is devoted to conclusions and to sketch future lines of research. 
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2 Circulant matrices and other related structures 



Let / be a trigonometric polynomial denned over the set Q = [0, 2tt) and having degree c > 0, i.e., 
f( u ) = Ylk=-c a k elku > i 2 = ~ !• From the Fourier coefficients of /, that is 



a q 



1 

2^ 



Q 



f(u)e- iju du, 



3 G 



(1) 



one can build the circulant matrix C n (f) = [o( r _ g ) modn]™ s l - For example, let /(u) = 3 — 2e m + 
e~ 2lu . The degree of / is c = 2 and we have ao = 3, a± = —2 and a_2 = 1; if we take n > (2c — 1) 



see the discussion below, since a 



(—2) modra 



03 the circulant matrix C${f) is given by 
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It is clear that the Fourier coefficient aj equals zero if the condition \j\ < c is violated. The 
matrix C n (f) is said to be the circulant matrix of order n generated by /, and can be written as 
C n (f) = X)m< c a jZh, where the matrix 




1 





1 


1 



is the cyclic permutation Toeplitz matrix. In addition, if F n denotes the Fourier matrix of size n, 
i.e. 

■ n-l 
j,k=0 ' 



1 



2irijk 

e ^~ 



(2) 



then it is well known (see e.g. [S]) that 

C n (f) = F n A n (f)F^ 

where 



j=l,...,n V J 
diag( v / nF ? ^a 



27TJ 

n 

[ao, ai, 



(3) 



(4) 



j "n-l 



a being the first column of the matrix C n {f). 

Under the assumption that c < \_(n — l)/2j , the matrix C n {f) is the Strang or natural circulant 
preconditioner of the corresponding Toeplitz matrix T n (f) = [o(r-«)] r s _ ( see E] an d the references 
therein). We observe that the above-mentioned assumption c < [_(rz — 1) /2j is fulfilled at least 
definitely, since each c is a fixed constant and n is the matrix order: in actuality, in real applications 
it is natural to suppose that n is large if we assume that C n {f) comes from an approximation process 
of an infinite-dimensional problem. However if the symbol / has a zero at zero (this happens in the 
case of approximation of differential operators), then C n (f) is singular and it is usually replaced by 
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We end this section with the definition of g-circulant matrix. A matrix C n;5 



a rank-one correction that forces invertibility: the latter is called in the relevant literature modified 
Strang preconditioner. 

of size n is called 

'n,g — [ a (r~gs) modn]" s _ : f° r an introduction and for the 
algebraic properties of such matrices refer to Section 5.1 of the classical book by Davis [8], while 
new additional results can be found in \18\ 126] and references therein. For instance if n = 5 and 
g = 3 then we have 



g-circulant if its entries obey the rule C n 



C n . 



a ci2 cla a\ a 3 

a\ a 3 ao a 2 04 

a 2 04 a\ a 3 ao 

a 3 ao a2 04 a± 

a 4 ai a 3 ao a 2 



Also in this case, as in ordinary circulant setting, the coefficients aj, j G Z, could arise from a 



given symbol / (see ([!])). For instance, with f{u) 
a_2 = 03 = 1 so that 



3 - 2e m + e~ 2tu , we find a = 3, m = -2, and 
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3 Two-grid and Multigrid methods 

Let A n G C raxn , and x n , b n G C n . Let G C nxfc , A; < n, a given full-rank matrix and let us 
consider a class of iterative methods of the form 



(5) 



where A n = W n - N n , W n nonsingular matrix, V n := I n - W' l A n G C nxn , and b n := W~ l b n G C 
A Two-Grid Method (TGM) is defined by the following algorithm: 

TGMfc, Cpost^n)^ ) 

0. X n = Vnfprc(^n \ &n,prc) _ _ " 

1. d n — A n x n b n 

2. 4 = (p^d n 

3. ,4 fc = 

4. Solve = dk 

5. x n — x Ti 

6. = 



HifpostO^U &n,post) 



Steps 1. — > 5. define the "coarse grid correction" that depends on the projector operator p n , 
while Step 0. and Step 6. consist, respectively, in applying u pre times and ^ pos t times a "pre- 
smoothing iteration" and a "post-smoothing iteration" of the generic form given in ([5]). The global 
iteration matrix of the TGM is then given by 



TGM(T4i l p re , ^i,post)Pn) — Ki,post 



In ~ Pn ((PT A nPn y\p k n ) H An 



*n,pre- 
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If k is large, the numerical solution to the linear system at the Step 4. could be computa- 
tionally expensive. In such case a multigrid procedure is adopted. Fix < m < n, the se- 
quence < n m < n m _i < ■ ■ ■ < rii < uq = n and the full-rank matrices pl^i £ C n '" lXn ', 

for i = l,...,m. The multigrid method produces the sequence {x„ }keN defined by Xn = 
MGM(K^,V^ p p ° o s s t t,K 1 ,^n,6n,^,0)(x^ ) ) with the function MGM defined recursively as follows: 

= MGM(Ce,0'^ +1 .\A^.i)(4 ) ) 

If i = m then Solve A ni Xn^ = b ni 
Else 

0- X ni = Vji^, pre^n, j ^rii,pre) 

1- d ni — A ni x ni b ni 

2. d ni+1 = { P n nn H d ni 

for s = 1 to 9 

Xn i+ i = yiGM.(y n ^ +1 ,pTe, Ki 4 +i,post)P«i+n ^"i+i i dn i+1 , 9, 1 + l)(£ni ^) 
^n-i — x ni Pru ^n i+1 
6- x rti = 1 / ' ra P p 0st (Xn i , ftni.post) 

The choices 9 = 1 and 9 = 2 correspond to the well-known V-cycle and W-cycle, respectively. 

In the present paper, we are interested in proposing such kind of techniques in the case where A n 
is a Toeplitz matrix. However, for a theoretical analysis, we consider circulant matrices according 
to the local Fourier analysis for classical multigrid methods (see [9]). Even if we treat in detail 
the circulant case, in the spirit of the paper [2], the same ideas can be plainly translated to other 
matrix algebras associated to (fast) trigonometric transforms. First we recall some convergence 
results from the theory of the algebraic multigrid method given in |19j . 

By || • ||2 we denote the Euclidean norm on C™ and the associated induced matrix norm over 
C nxn . If X is positive definite, || • ||x = ||^ 1//2 • H2 denotes the Euclidean norm weighted by X on 
C n and the associated induced matrix norm. Finally, if X and Y are Hermitian matrices, then the 
notation X < Y means that Y — X is nonnegative definite. In the following we use some functional 
norms: more precisely the usual L°° norm || • ||oo defined 

II /II 00 — s^Pxgo 1/(2;) I, and the weighted 
L 1 norm || • ||i defined as ||/||i = ^ Jq \f{x)\dx (according to the Haar measure). 

Theorem 3.1 ([E]). Let A n be a positive definite matrix of size n and let V n be defined as in the 
TGM algorithm. Suppose that there exists a pos t > independent of n such that 

\\Vn ^ost^n \\a„ — ll^^llyln C^post \\%n || A n D~ 1 A n ' 

Vx„ G C n , (6) 

where D n is the main diagonal of A n . Assume that there exists 7 > independent of n such that 

min \\x n -PnH\\ 2 D n < 7lWli„, Vx n G C". (7) 



y<EC k 

Then 7 > a post and 



1 1 TGM (I, ^pos t , Pn) I Un < A/ 1 - "pWt- 
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Conditions ([6]) and (J7]) are usually called as "smoothing property" and "approximation prop- 
erty", respectively. 

We note that a post and 7 are independent of n and hence, if the assumptions of Theorem 13.11 
are satisfied, the resulting TGM is not only convergent but also optimal. In other words, the 
number of iterations in order to reach a given accuracy e can be bounded from above by a constant 
independent of n (possibly depending on the parameter e). 

Of course, if the given method is complemented with a convergent pre-smoother, then by the 
same theorem we get a faster convergence. In fact, it is known that for square matrices A and B 
the spectra of AB and BA coincide. 

Therefore TGM(y^p r c c , l^p^ t ,pj) and TGM(J, V^p r c e ^p° S g t ,p^) have the same eigenvalues so 
that 

||TGM(V^pre, ^post)Pn)IUn = l|TGM(J,y^e< P p ^ t ,P^)|U n < ^ " «post/7, 

and hence the presence of a pre-smoother can only improve the convergence. 

Concerning multigrid methods, in [19] the V-cycle convergence is considered with a result which 
could be seen as the analog of Theorem 13.11 For other bounds concerning the convergence rate of 
the V-cycle see [T7] and reference therein. Regarding the convergence of the W-cycle, we point out 
that a rigorous TGM analysis is sufficient for determining the optimality of the W-cycle (see [27J). 



4 Projector operators for circulant matrices 

Let A n := C n (f) be a circulant matrix generated by a trigonometric polynomial /. In order to 
provide a general method for obtaining a projector operator from an arbitrary banded circulant 
matrix P n , for some bandwidth independent of n, we introduce the operator Z^ g G R nxfc , k = ^ G 
N, where 



1 if r = (mod n) 
otherwise, 



i = 0, . . . ,n — 1, 
j = 0, . . . , k - 1. 



(8) 



n,g 

The operator Z^ g represents a special link between the space of the frequencies of size n and the 
corresponding space of frequencies of size k. 

Lemma 4.1. Let F n be the Fourier matrix of size n defined in (TJ) and let Z^ g G R nxk be the 

(9) 



matrix defined in (0) • If k = ~ G N then 



F^Zl 



1 



n n,g 



V9 n ' 9 k ' 



where I n „ G 



nnxk 



and 



L n,g 



\ h ' 






h 












[ Ik _ 







g times, 



with Ik being the identity matrix of size k. 
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This simple relation (see \18\ Lemma 3.3 and Remark 3.5] for the details of the proof and [26] 
for recent findings on these structures) is the key step in defining an algebraic multigrid method, 
since it allows us to obtain again a circulant matrix at the lower level. Indeed, denoting by A n the 
diagonal matrix obtained from the eigenvalues of A n (see flU)), we infer that := -i^ g A n J n)9 is 
again a diagonal matrix. Therefore 

7-1 jT A J T?H 

— g r k 1 n,g L -*n 1 n,g r k 

= A k , 

where A k is a new circulant matrix. Consequently, starting from the matrix Z!* g , it is possible to 
define a generic projector 

Pn,g = PnZ n ,gi (10) 

where P n is a circulant matrix. Indeed A n P n is a circulant matrix and then A k = (Pn, g ) H ^nPn,g 
is again a circulant matrix. We note that, since k = G N, n must be a multiple of g. We are left 
to determine the conditions to be satisfied by P n = C n (p) (or better by its generating function p) 
in order to get a projector which is effective in terms of convergence. 

Definition 4.1. Given x G [0, 27r) ; g 6 N, g > 2, the set of ^-corners of x is £l g (x) = {y = x + 
(mod 2ir) \ k = 0, . . . , g — 1} and the set of g-mirror points is Ai g (x) = £l g (x) \ {x}. 

TGM conditions Let A n := C n (f) with f nonnegative, trigonometric polynomial and let p\ g = 
C n (p)Z^ g with p trigonometric polynomial. Assume that f(xo) = for xq G [0, 2tt), choose p such 
that the following relations 

hm^M<oo, \/yeM g (x), (11) 

x-+xo f(x) 

P 2 (V)>0, VxG[0,2vr), (12) 

yeVt g (x) 

are fulfilled. 

If / has a unique zero xq G [0, 2n), then we set P n = C n (p) where p is a trigonometric polynomial 
defined as 

p( x )= J] (2-2cos(x-x)) r/3/21 ~ Yl \x-x\ 2 ^ /2 \ (13) 

X£Mg(x ) X€Mg(x ) 

for x G [0, 27r), with 

2i 



P > Anin = min 



hm — — - — < +oo 

x->x f (x) 

thus conditions (jlip and (|12p are satisfied. 

Before proving (in Section I5.ip that conditions (jlip and (I12p are sufficient to assure the TGM 
optimality, we consider a crucial result both from a theoretical and a practical point of view. 
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Proposition 4.1. Let f be a nonnegative function, k = — G N, p\ Q = C n (p)Z k G <C nxk , with p 

_____ 5 '** 

trigonometric polynomial satisfying the condition (|77|) /or any zero o/ / and globally the condition 
(fi^j) . T/ien i/te matrix (Pn, g ) H C n (f)Pn,g C fcxfc coincides with Ck(f) where f is nonnegative and 



(14) 



2/en s ( 



/or x G [0, 27r), i.e., i/ie projected matrix is obtained picking every g-th entry out of the symbol 
f\p\ 2 . In particular 

1. if f is a polynomial then f is a polynomial with a fixed degree - , where c is the degree of 

f\v\ 2 ; 

2. if xq is a zero of f then f has a corresponding zero yo where yo = gxQ ( mod 2ir ); 

3. the order of the zero yo of f is exactly the same as the one of the zero xq of f, so that at the 
lower level the new projector can be easily defined in the same way. 

Proof. First we observe that, from ([3]) and (jJJ), 

{Pic,)" Cn{f) P i g = (Zl g ) H (C n {p)) H C n (f)C n {p)Zl g 

= (Z^) H C n (f\p\ 2 )Z^ g , 

thus the generating function of the circulant matrix (C n (p)) H C n (f)C n (p) is f\p\ 2 ■ Denoting by atj 
the Fourier coefficients of f\p\ 2 , it holds 



C n (f\p\' 



m-l 



(l — s) modnJ rs _g> 



and then, by ([8]), the entries of the matrix (Z% „) H C n (f\p\ 2 )Z k are given by 



Cn(/M 2 )< s 



n-1 



zZ{Cn(f\p\ 2 )]rAZn, : 



1=0 
n-1 

E< 



''(r—() modn' 



(a) ^(r—gs) modni 



r = 0, . . . , n — 1, S = 0, 



k-h 



(z k n j H c n (f\p\ 2 )z*, : 



n-1 

Y^ Z i 9 ) H U\Cn{f\p\ 2 )Z h n, 9 ]^ 

e=o 

n-1 

^ ] &£—gr a (£—gs) modn 



(b) ^(gr—gs) modn) 



r, s = 0, . . . , k — 1, 



where (a) follows because there exists a unique £ G 1,2,... , n — 1 such that £ — gs = (mod n), 
that is, £ = gs (mod n) and, since < gs < n — 1, we obtain £ = gs; similarly for (b). Now 
if we denote by Oj the Fourier coefficients of / it only remains to show that 6( r _ c ) = a^ gr _ gc ^, 
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r, c = 0, . . . , k — 1, from which we directly infer that [C n (/)] rjC = 6( r _ c ) mo dn = a (gr~gc) modn- Since 
/|p| 2 is a polynomial, we can always write 

oo oo 

f\p\ 2 (x) = E a ^> />) = E ( 15 ) 

£=— oo £=— oo 

Prom ©, ((TU) and ([15]) . we have 

- 70 y j=0£=-oo 



Recalling that 



lv-i «*a \ \ \ii = gt , 1 f 2w ilxj f I if r - () 

- > e 3 = S n ^ and — / e tlx dx = < n 



a ^-^ 1 otherwise 2tt n 10 otherwise ' 



we find that 



+oo „ 2 7r 

E I j x{t - {r - c)) dx 



t=-oo — J0 

= a g(r-c)- (16) 

From the expression of /, since f(xo) = then it must hold p(y) = Vy G 7W 9 (xo) to satisfy 
the (jlip . Thus yo = y^o (mod 2ir) is a zero of / (i.e. item 2. is proved). 

Moreover, by (fT2|) . we deduce that j9 2 (xo) > since p 2 (y) = 0, Vy € .A/f 9 (xo), and the order 
\2 f x 



of the zero yo of / l~j \p\ is the same as the order of f(x) at xq. Furthermore, by (jlip we 

can see that \p\ 2 ( x+ ^ k ^ has at yo a zero of order at least equal to the one of f(x) at xq, for any 

k = 1, . . . , g — 1. Since all the contributions in / are nonnegative the thesis of item 3. follows. 

Finally we have to prove item 1. Since bj are the Fourier coefficients of / and cij are the Fourier 
coefficients of the polynomial f\p\ 2 , see (fT5j) . from (fT6j) we deduce that 

ev" E "<"'""• 



□ □ 



Hence, if the polynomial has degree c, / has degree at most 

5 Proof of convergence 

Using the results in Section U it is possible to prove the optimality of the TGM and of the W-cycle 
(W-cycle requires g > 2). 
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5.1 TGM convergence 

The smoothing property for g = 2 was proved in [23] and it holds unchanged also for g > 2. 

Lemma 5.1 (|23|). Let A n := C n (f) with f being a nonnegative trigonometric polynomial (not 
identically zero) and let V n := I n — ujA n , 0<w<2/||/|| oo . If we choose a pos t so that a pos t < 
aow(2 — wH/Hoo) then relation © holds true. 

If in the previous Lemma we choose ui = \\f\\~^, then a pos t < ||/||i/||/||oo and the best value of 
a post is a pos t,best = ||/||i/||/||oo- Moreover, the result of Lemma [5TT1 can be easily generalized when 
considering both pre-smoothing and post-smoothing as in [TJ. 

The following result shows that TGM conditions (jlip and (|12p are sufficient in order to satisfy 
the approximation property. 

Theorem 5.1. Let A n := C n (f) with f being a nonnegative trigonometric polynomial (not iden- 
tically zero) and let p k g = C n {p)Z k g the projector operator, with Z k g defined in and with p 
trigonometric polynomial satisfying condition ([77]) for any zero of f and globally the condition (| 
Then, there exists a positive value 7 independent of n such that inequality Q holds true. 

Proof. The proof is similar to that of Lemma 8.2 in [22], but we report it here for completeness. 
First, we recall that the main diagonal of A n is given by D n = aoI n with oq = (2ir)~ 1 Jq f = 
H/lli > 0, so that || • \\ 2 Dn = a || • |||. 

In order to prove that there exists 7 > independent of n such that for any x n G C" 

min \\x n -p k yf Dn = a min \\x n -p k n y\\ 2 2 < j\\x n \\ 2 A , 

we chose a special instance of y in such a way that the previous inequality is reduced to a matrix 
inequality in the sense of the partial ordering of the real space of the Hermitian matrices. For any 
x n 6 C n , let y = y(x n ) 6 C k be defined as 



!J 

Therefore, (|7|) is implied by 



/ k \H k 

\Pn,g) Pn,g 



-1 



(Pn,n) x n- 



\\ x n-Pn, g y\\l < (7/«o)lkn||i n , Vx„ E C n , 
where the latter is equivalent to the matrix inequality 

W n (p) H W n (p) < ( 7 /a )C n (/), (17) 

with W n {p) = I — Pn g [(Pn, g ) H Pn,g] 1 (Pn,g) H ' • Since, by construction, W n {p) is an Hermitian 
unitary projector, it holds that W n (p) H W n (p) = W 2 (p) = W n (p). As a consequence, inequality 
(1171) can be rewritten as 



to 



W n {p) < (7/a )C n (f). (18) 
, ,. ^ composition in ©, p^ g 



If k = ^ € N, following the decomposition in Q, p\ „ = C n (p)Z k can be expressed according 



« 9 ) H = ^[A(°)| A «|...|Ar i) 



F H 
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where 



AW 

p 



0,...,g-l, 



» 



with x y - L> = 2-nj /n. 

Let £ C 9 whose entries are given by the evaluations of p over the points of f2(x^ 
for = 0, . . . , k — 1. There exists a suitable permutation by rows and columns of F^W n (p)F n , 
such that we can obtain a g x g block diagonal matrix whose ^th diagonal block is given by 
I g — p[[A(p[[j]) T /WpifAM- Therefore, using the same notation for f[fi] and denoting by diag(/[/x]) 
the diagonal matrix having the vector f[n] on the main diagonal, the condition (|18|) is equivalent 
to 

p[fA(p[fA) T 



In 



\\pM\1 



< (7/oo)diag(/Dz]), 



(19) 



for /j, = 0, . 
of 



, k — 1. By the Sylvester inertia law |12j . the relation (|19p is satisfied if every entries 



diag(/[ / x])- 1 /2 j 



pM(pM) rj 
IIpMIII 



diag (/[//]) 



-1/2 



is bounded in modulus by a constant, which follows from the TGM conditions (jlip and ([12 
Furthermore, if we put 

p 2 (y) 



max 



/(*) 



the condition ([7]) is satisfied choosing a value of 7 such that 7 > g(g — l)aohz. 



□ 



□ 



Combining Lemma f5.1l and Theorem l5.1l with Theorem l3.ll it follows that the TGM convergence 
speed does not depend on the size of the linear system. 



5.2 Multigrid convergence 

The optimal TGM convergence rate proved in Theorem 15.11 can be extended to a generic recursion 
level of the multigrid procedure obtaining the so called "level independency" property. The key 
tools to do that are the Proposition 14.11 and an explicit choice of the projector, considering for 
instance the symbol p in (|13p . Indeed, the "level independency" was already proved in literature 
for g = 2 (see [U HJ E]) and the proof can be extended to g > 2, as in Theorem 15.11 

The "level independency" implies that the W-cycle has a constant convergence rate independent 
of the problem size [27] . However, the fact that the convergence speed does not depend on the 
size of the linear system does not implies the optimality of the method, because the computational 
work at each iteration is not taken into account. For estimating the computational work at each 
iteration of a multigrid method, we have to consider the size of the coarse problem and the number 
8 of recursive calls. In our case the size of the problem at the level i is rii = grii-\. According to the 
analysis in |27j . we assume that the multigrid components (smoothing, projection, . . . ) require a 
number of arithmetic operations which is craj, with c constant independent of rij, up to lower order 
term. From equation (2.4.14) in [27], the total computational work C m of one complete multigrid 
cycle is 



Cm 



jcn 



0(n log n) 



for 9 < g 
for 6 = g 



(20) 
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where the symbol = means equality up to lower order terms. It follows that for g = 2 the W-cycle 
can not be optimal even in the presence of "level independency" , because each multigrid iteration 
requires a computational cost of 0(n log n) while the matrix vector product is of 0(n). On the 
other hand, for g = 3 the W-cycle has C m = 3cn and hence it is optimal if the "level independency" 
is satisfied. More in general, the proposed multigrid will be optimal for a number G N of recursive 
calls such that 1 < 9 < g. 



5.3 Some pathologies eliminated when using our algorithm 

From conditions (3.4) and (3.5) in [23], we know that, for g = 2, if xq is a zero of /, then /(xq + ir) 
must be positive: otherwise relationship (3.5) in [23] cannot be satisfied with any polynomial p. 
But if we consider g = 3 this is no longer a problem, because conditions (jlip and (|12p impose 
that, if xq is a zero of /, then / (rco + §7r) and / (xq + §7r) must be positive, while there are no 
conditions on f(xo + n). 

For g = 3, if / has a unique zero xq G [0, 2tt) of finite order, then we consider x = xq + |7r 
(mod 27r) and x = xo + |vr (mod 27r) and we set P n = C n (p) where p is a trigonometric polynomial 
defined as 

(21) 

p(x) = (2-2cos(x-x))^/ 2 l(2-2cos(x-5))^/ 2 l ~ \x - x\ 2 ^^\x - x\ 2 ^/ 2 \ 
for x G [0, 27r), with 

|2i 



P > Anin = min 



hm — - — < +oo 

x^xq j(x) 



thus conditions (jlip and (|12p are satisfied. If / shows more than one zero in [0, 27r), then we 
consider a polynomial p which is the product of the basic polynomials of kind (|2ip . satisfying the 
condition (jlip for any single zero and globally the condition f)12|) . 

Example: The symbol 

f(x) = (2-2cos(x))(2 + 2cos(x)) 

vanishes at zero and tt with order two. For g = 3, we have A^3(0) = {^, ^} and Ai^i^) = 
{if' T"}' ^ nus the trigonometric polynomial 

p(x)= Yl (2-2cos(x-x)) (22) 

x£M 3 (0)\JM 3 (n) 



satisfies the TGM conditions (|lip and (I12p and defines an optimal TGM. 



6 Numerical experiments 

In this section, we apply the proposed multigrid method to symmetric positive definite circulant 
and Toeplitz systems A n x = b. We choose as solution the vector x such that X{ = i/n, i = 1, . . . , n. 
The right-hand side vector b is obtained accordingly. As smoother, we use Richardson with cjj = 
l/ll/jlloo) f° r j = 0, . . . , m — 1 (m is number of subgrids in the algorithm, m = 1 for the TGM), for 
pre-smoother and the conjugate gradient for post-smoother. In the V-cycle and W-cycle procedure 
when the coarse grid size is less than or equal to 27, we solve the coarse grid system exactly. The 
zero vector is used as the initial guess and the stopping criterion is 1 1 t™^ 1 1 2 / 1 1 1 1 2 — W 7 , where r q 
is the residual vector after q iterations and 10 is the given tolerance. 
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6.1 Cutting operators for Toeplitz matrices 

When dealing with circulant matrices, using the projector defined in (|1Q|) . the matrix at the lower 
level is still a circulant matrix, while for Toeplitz matrices, if we consider A n := T n (f) and p\ 3 = 
T n (p)Z^ 3 , where p is defined in accordance with the formula (f2T|) and fc = | S M, we find that 

T n (p)T n (f)T n (p) = T n (fp 2 ) + G n (f,p). 

Furthermore, if 2/3 + 1 is the bandwidth of T n (p), the matrix G n (f,p) has rank 2/3 and is formed 
by a matrix of rank (5 in the upper left corner and a matrix of the same rank in the bottom right 
corner. According to the proposal in [2], we take a cutting matrix that will completely erase the 
contribution of G n (f,p), so that, at the lower level, the restriction of the matrix T n (p)T n (f)T n (p) 
is still a Toeplitz matrix and thus we can recursively apply the algorithm. The proposed cutting 
matrix is as follows: 

Z~k 
n,3 — 

where 0^ _r is the zero matrix of size j3 x (k — r); Z k 3 has the first and the last (3 rows equal to zero 
and therefore it is able to remove corrections of rank less than or equal to 2/3. Since G n (f,p) has 
rank 2/3, we deduce that A k _ r = (Z% 3 ) H T n (p)T n (f)T n (p)Z* 3 is Toeplitz and we cannot obtain a 
Toeplitz matrix of size greater than this. As a consequence, for Toeplitz matrices, the projector is 
then defined as 

p£, 3 = r„(p)z* 3 . 

Also the size of the problem should be chosen in such a way that a recursive application of the 
algorithm is possible; in our case, if we choose n = 3 a — £ with £ = /3 — 1, then the size of the problem 
at the lower level becomes k! = k - r = n ~ 2( 3 /3 ~ 1) = 3 "~ (/3 ~ 1 ; ]~ 2(/3 ~ 1) = - (fi - 1) = - £. 

6.2 Zero at the origin and at n. 

We present some examples where the generating functions /o vanish at the origin and at ir. Firstly, 
we consider the Example 15.31 where the symbol 

/ (x) = (2 - 2cos(»)(2 + 2cos(»), 

vanishes at and it with order 2. According to (f2T|) . we choose the projector p^ 3 = C n (po)Z^ 3 if 

A n is a circulant matrix and p k n 3 = T n (po)Z k 3 if A n is a Toeplitz matrix, where Po = P defined in 

(|22p . Fixing = and = it, the position of the new zeros \ for k = 1, 2, . . . , m — 1 with 
j = 1,2, move according to Proposition 14.11 and, in this case, the functions p k are equal to p for 
every level k. Tables [Hand [2] report the number of iterations required for convergence in the case of 
circulant and Toeplitz systems, respectively. In all cases we note an optimal behavior at exception 
of the V-cycle for Toeplitz matrices where the number of iterations slightly grows with the size n. 
In the second example we increase the order of the zero in it considering the function 

f (x) = (2 - 2cos(x))(2 + 2cos(x)) 2 . 

which has a zero at with order 2 and one at n with order 4. The polynomial po = P defined in 
(|22p still satisfies the TGM conditions (jlip and f)12|) . The functions p^ do not change at the lower 



Z 



2(0-1) 



n-2/3,3 
k—r 



ri Y k — r 
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Table 1: Circulant case: fo(x) = (2 — 2cos(x))(2 + 2cos(x)). 



n 








# iterations 










Two- 


grid 


V-cycle 




W-cycle 






^pre — 


^prc — 


^pre — ^pre — 


^pre 


k'pre — 






^post = 1 


^post = 2 


^post — 1 ^post — 2 


^post : 


= 1 ^post = 2 


3 4 


= 81 


11 


6 


11 6 


11 


6 


3 5 


= 243 


11 


6 


11 7 


11 


6 


3 6 


= 729 


11 


6 


11 7 


11 


6 


3 7 


= 2187 


11 


6 


11 7 


11 


6 



Table 2: Toeplitz case. / (x) = (2 - 2cos(x))(2 + 2cos(x)) 



n 










# iterations 










Two- grid 


V-cycle 


W-cycle 








^pre — 


^pre — 


^pre = ^pre = 


^pre = ^pre 








Impost = 1 


^post = 2 


^post = 1 ^post = 2 


: = 1 ^post = 2 


3 4 


- 3 = 


78 


24 


14 


24 14 


24 14 


3 5 


-3 = 


240 


24 


15 


35 20 


28 16 


3 6 


-3 = 


726 


24 


15 


43 24 


29 16 


3 7 


-3 = 


2184 


24 


15 


49 27 


29 16 
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Table 3: Circulant case. / (a?) = (2 - 2cos(x))(2 + 2cos(a;)) 2 , tolerance=l(r 3 



n 








# iterations 








Two-, 


grid 


V-cycle 


W-cycle 






^pre — 


^pre — 


^pre — ^pre — 


^pre — ^pre — 






t'post = 1 


^post = 2 


^post — 1 ^post — 2 


^post — 1 ^post — 2 


3 4 


= 81 


20 


9 


20 9 


20 9 


3 5 


= 243 


20 


9 


18 9 


20 9 


3 6 


= 729 


20 


9 


18 9 


20 9 


3 7 


= 2187 


20 


9 


18 9 


20 9 



Table 4: Toeplitz case. f (x) = (2 - 2cos(x))(2 + 2cos(x)) 2 , tolerance=10~ 3 



n 










# iterations 










Two- 


grid 


V-cycle 


W-cycle 








^pre — 


^pre — 


^pre — ^pre — 


^pre — ^pre — 








^post = 1 


^post = 2 


^post — 1 ^post — 2 


^post — 1 ^post — 2 


3 4 


- 3 = 


78 


50 


31 


50 31 


50 31 


3 5 


- 3 = 


240 


48 


31 


93 35 


72 32 


3 6 


- 3 = 


726 


47 


31 


74 34 


68 31 


3 7 


- 3 = 


2184 


47 


31 


76 34 


68 31 



levels. In Tables [3] and H] we report the number of iterations required for convergence in the case of 
circulant and Toeplitz systems, respectively. Since /o has a zero of order 4 the condition number 
of A 3 7 = 0((3 7 ) 4 ) = O(10 13 ). Therefore, using double precision, for this example we choose a 
tolerance equal to 10~ 3 . This choice agrees also with the plots in Figures [T] and [2] where we note 
an optimal reduction of the residual norm only until about 10~ 3 . 

The last example of this subsection is taken from [5]. The generating function 

fo(x) = 6 - 4cos(2x) - 2cos(4x), 

vanishes at and it with order 2. The symbol of the projector is again po = P defined in (|22p . 
The initial guess is a random vector u such that < uj < 1, the pre-smoother is a step of damped 
Jacobi with parameter uij = [A/J^i/H/^lloo while the post-smoother is a step of damped Jacobi 
with parameter ujj = 2[A ? -]i i i/||/(x)|| 00 for j = 0, ... ,m — 1. The coarser problem is fixed such 
that is has size lower than 6. Table [5] shows that the number of iterations required to achieve the 
tolerance 10 -7 remains constant increasing the size n of the system like for the multigrid technique 
proposed in [5]. The number of iterations is reasonable in both cases even if a direct comparison 
can not be done because of the difference in the choice of the projection techniques and in the size 
of the projected problems. 
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Figure 1: Circulant: Graph of the residual in logarithmic scale of the V-cycle (left) and W-cycle 
(right) with different sizes n, with f pre = z^ post = 1 and a fixed number of iterations iter = 400; 
fo(x) = (2 - 2cos(x))(2 + 2cos(x)) 2 . 




Figure 2: Toeplitz: Graph of the residual in logarithmic scale of the V-cycle (left) and W-cycle 
(right) with different sizes n, with f pre = z^ post = 1 and a fixed number of iterations iter = 400; 
f (x) = (2 - 2cos(x))(2 + 2cos(x)) 2 . 



Table 5: Toeplitz case. fo(x) = 6 — 4cos(2x) — 2cos(4x), = u post = 1, tolerance=10 



n 








# iterations 










Two-g 


rid W-cycle 


V-cycle 


3 4 


-3 = 


78 


15 


19 


28 


3 5 


-3 = 


240 


15 


20 


39 


3 6 


-3 = 


726 


14 


20 


45 


3 7 


-3 = 


2184 


13 


20 


47 
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Table 6: Toeplitz case. fo(x) = (2 — 2 cos [x — tolerance=10 



11 








# iterations 










V-cycle 


W-cycle 








^pre 


^pre — 


^pre — ^pre — 








^post = 


= 1 ^post = 2 


^post — 1 ^post — 2 


3 4 


- 1 


= 80 


33 


37 


33 37 


3 5 


- 1 


= 242 


30 


31 


30 31 


3 6 


- 1 


= 728 


30 


31 


30 31 


3 7 


- 1 


= 2186 


30 


31 


30 31 



6.3 Some Toeplitz examples 

In this subsection we consider only the more interesting case for practical applications: Toeplitz 
matrices with a multigrid strategy. 

The first example is a function with a zero not at the origin or ir: 

fo(x) = (2 -2 cos (x- I 

which vanishes at ir/3 with order 2. Moreover, we choose as true solution a random vector instead 
of a smooth solution. The tolerance is again 10 -7 . The symbol of the projector at the first level is 

Po{x) = (2 — 2 cos (x — 7r)) ^2 — 2 cos ^x — -ir 

while at the lower levels it changes with the zero of fj which moves according to Proposition 14.11 
Table [6] shows an optimal convergence both for V-cycle and W-cycle. 

In the second example, we consider the dense Toeplitz matrix generated by the function f(x) = 
x 2 , which has the Fourier series expansion 

„, . 7r 2 / cos(x) cos(2x) cos(3x) 
/(x) = --4 — ^-—^ + 



3 \ l 2 2 2 3 2 

Such function shows a unique zero at with order 2 and hence we use the projector with symbol 

Po(x) = ^2 — 2 cos ^x — -iT^j ^ ^2 — 2 cos ^x — -ir 

In Table [7] we report the number of iterations required for the convergence with the preassigned 
accuracy and we note again the optimal behavior. 

7 Conclusions and future work 

In this paper we have extended the rigorous two-grid analysis for circulant matrices to the case 
where the size reduction is performed by a factor g with g > 2. The interesting novelty is that 
the new size reduction strategy allows to eliminate some pathologies which occur when g = 2. In 
particular, if the considered matrices come from the approximation of certain integro-differential 
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Table 7: Toeplitz case. f(x) = x 2 



n 








# iterations 










V-cycle 


W-cycle 








^pre 


^pre — 


^pre — ^pre — 








^post = 


= 1 ^post = 2 


^post — 1 ^post — 2 


3 4 


- 1 


= 80 


21 


11 


21 11 


3 5 


- 1 


= 242 


18 


11 


21 11 


3 6 


- 1 


= 728 


18 


11 


21 11 


3 7 


- 1 


= 2186 


18 


11 


21 11 



equations then we have two source of ill-conditioning and the zeros of the underlying symbol are 
located at zero and at 7r: this situation is a special case of mirror point zeros and, when g = 2, 
it is possible to prove that the resulting two- grid iteration cannot be optimal (see [101 122|). Such 
difficulty can be overcome when we choose a larger g. Moreover, when increasing g the size of the 
coarse problems decreases: as a consequence more multigrid recursive calls could be considered, 
like the W-cycle which is proved to be optimal for g > 3. 

We stress that the numerical experiments are encouraging not only for circulant matrices but 
also regarding Toeplitz matrices and concerning the V-cycle algorithm. A future line of research 
must include the multilevel setting, following the approach in [221 [I], and a rigorous proof of 
convergence for the whole V-cycle procedure in accordance with the proof technique introduced in 

®- 
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